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Computation of Neutron Star Surface Emission Spectra 
for Arbitrary Magnetic Field Directions without Diffusion 
Approximation 
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Abstract. To derive physical properties of the neutron 
star surface with observed spectra, a realistic model spec- 
trum of neutron star surface emission is essential. Lim- 
ited by computing resources, a full computation of the 
radiative transfer equations without the diffusion approx- 
imation has been conducted up to date only for the case 
of local magnetic fields being perpendicular to the stel- 
lar surface. In this paper we report the full-computation 
result for an arbitrary field direction. For comparison we 
also compute the radiative transfer equation using the dif- 
fusion approximation. For a given effective temperature, 
the computed spectrum with the diffusion approximation 
is always softer than that of a full computation at a non- 
negligible level. It leads to an over-estimate of the effective 
temperature if the diffusion approximation spectrum is 
employed in the spectral fitting. Other characteristics for 
different magnetic field orientations, such as the beaming 
pattern of the two polarization modes and the structure 
of the atmosphere, are also discussed. 



1. Introduction 

The spectra of some radio pulsars and radio-quiet isolated 
neutrons reveal a thermal origin in the soft X-ray band 
(e.g. Becker & Triimper 1997; Pavlov, Zavlin & Sanwal 
2002). These thermal emissions are believed to originate 
from the surface of neutron stars and they carry the infor- 
mation of physical properties of the neutron stars. Mod- 
eling the atmosphere of neutron stars to obtain a more 
realistic spectrum of the surface emission to compare with 
observations may help us to have better knowledge about 
those properties such as the composition, the mass-to- 
radius ratio and the strength of the magnetic field at the 
surface. Many neutron star atmosphere models have been 
developed in various physical conditions (see e.g. Zavlin & 
Pavlov 2002, for a review). We construct an atmosphere 
model with pure hydrogen gas. In the field range of 10 11 - 
10 13 G and the effective temperature of several million 
Kelvin, the hydrogen gas can be well approximated as a 
fully ionized plasma (Lai & Salpeter 1997). Therefore only 



the scattering and the free-free processes are involved in 
our computation. We solve the radiative transfer equation 
of angle-dependent specific intensity for the two polariza- 
tion modes in a plasma. The magnetic field direction is 
allowed to be arbitrary. Up to date, such a full compu- 
tation without the diffusion approximation has been done 
only for the case of the magnetic field being perpendicular 
to the stellar surface. In the following sections we briefly 
describe the method we use and then discuss the results 
for different field orientation and the difference between 
the full computation and that with the diffusion approxi- 
mation. 

2. Assumptions and the radiative transfer 
equations 

The equation of state adopted for the fully ionized hy- 
drogen plasma is that of an ideal gas, P ~ 2p , where 
P, p, and T are the pressure, density and temperature 
of the hydrogen gas respectively; m p is the proton mass. 
Because the scale height of the neutron star atmosphere 
(~ 1cm) is much smaller than the radius of the neutron 
star (~ 10 6 cm), the atmosphere is treated as a plane- 
parallel surface. We assume the atmosphere is in hydro- 
static equilibrium, and all physical quantities arc indepen- 
dent of time. The hydrostatic equilibrium in the atmo- 
sphere is described by the Oppenheimer-Volkoff equation, 
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where z is the height in the atmosphere, M is the mass of 
the neutron star, G is the gravitational constant, and R is 
the stellar radius. The terms in the first two parentheses 
on the right hand side of equation (1) are very close to 
unity. Equation (1) can be written as 
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Thomson scattering opacity. 
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The transport of photons is described by the radiative 
transfer equation, which is 

jjg _ K SC + K jj ji _ K ff By 

f dr 

where the specific intensity (#r, 4>r) is a function of di- 
rection specified by the angle #r between the surface nor- 
mal and the propagation of the photon and the azimuthal 
angle </>r (</>r = 0° in the direction of the magnetic field). 
The two polarization eigen modes are denoted with the 
superscript i with i = 1 for the extraordinary (X) mode 
and i = 2 for the ordinary (O) mode. The other notations 
are that /i = cos 9r, t is the Thomson optical depth, B v is 
the Planck function, k 1 sc is the Thomson scattering opac- 
ity for the i-mode photons, and rejj = Kg(l — exp (-rgpr)) 
is the reduced opacity of the thermal bremsstrahlung for 
the i-mode photons. All the opacities can be found, for ex- 
ample, in Meszaros (1992). The two boundary conditions 
for the above equation are 

r v {T = 0,fi,^ R ) =0 for - 1 < a* < 0, 
1 <\B 

-C( r = OO,^,0 R ) = -{B v +M-T -)|r=oo 

for < (j, < 1. (4) 

In our computation, close to the lower boundary deep in 
the atmosphere the frequency of low-energy photons can 
be lower than the local plasma frequency. In such a case, 
the O-mode photon has a non-zero imaginary part in the 
refraction index. To deal with such a situation, we define 
an additional opacity ks to represent the opaqueness of the 
plasma with ks = 2NlUJ , where Ni is the imaginary part of 

pc 

the refraction index N of the O-modc photons (Rajagopal 
et al. 1997). In the radiative equilibrium assumption, the 
total radiative flux is constant, that is, 4^ = 0, where 
H = Ef=i/f4M^/(47r). 
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Fig. 1. Temperature profiles of model atmospheres with 
g* = 10 14 cm/s 2 , T c = 10 6 K, and B = 10 12 G (the solid 
line is for the case without magnetic fields). 6b is the angle 
between the magnetic field and the surface normal. 



two boundary conditions. We choose the gray-atmosphere 
temperature profile, T 4 = |T 4 (r + |), as the first trial 
profile to start the iteration to obtain the radiative trans- 
fer equation solution in radiative equilibrium for a given 
effective temperature. In our computation we choose 20 
logarithmically equally-spaced grid points in a decade of 
the Thomson optical depth for the range from r = 10~ 5 
to t = 10 3 and similarly 10 grid points per decade in 
the frequency space for the photon frequency from 10 16 to 
10 19 Hz. Typically specific intensities in about 90 different 
directions in a hemisphere are involved in the computa- 
tion for an arbitrary magnetic field orientation. The num- 
ber of directions is significantly reduced when the field is 
perpendicular to the stellar surface because of the axial 
symmetry. The temperature correction is performed with 
the Unsold-Lucy procedure (Mihalas 1978). In a typical 
computation, about 100 iterations are needed to achieve 
an accuracy at the level of 1%. 



3. Numerical computation 

We use the Fcautrier method (Mihalas 1978) to solve the 
radiative transfer equation and also adopt the improved 
Feautrier method for reducing numerical round-off errors 
(Rybicki & Hummer 1991). The mean intensity-like and 
the flux-like quantities are defined as 

p i {k) = Ij(k) +Il(-k) ^ 

K(k) = Il{h) ^ ( ~ fc) - (5) 

where k is the unit wave vector of the out-going photons. 
With the above formulation, the radiative transfer equa- 
tion becomes a second-order differential equation with 



The diffusion approximation 

We also conduct computations with the diffusion approx- 
imation (e.g. Shibanov et al. 1992; Ho & Lai 2001), in 
which the specific intensity is approximated as P v = 
j^{u\, + -liFl) at all depth, where u % v and are the en- 
ergy density and flux density for the z-mode respectively. 
With this relation the radiative transfer equation for the 
specific intensity can be turned into that for the energy 
density, which can be solved with much less computing 
resource. However, it is obvious that in the presence of 
magnetic fields, the specific intensity cannot be well ap- 
proximated by this form (see Fig. 2). The same as in the 
case of the full computation, we use the Feautrier method 
to solve the radiative transfer equation and the Unsold 
Lucy procedure for the temperature correction iteration. 
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Fig. 2. The beaming pattern as a function of the photon 
propagation polar angle 0r. This figure shows the beaming 
pattern of a 10 17 -Hz photon for different field orientation 
and polarization modes. Parameters used are the same as 
those in Fig. 1 and </>r = 0°. 
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Fig. 3. Surface emission intensity spectra obtained with 
the full computation for different #r and different polar- 
ization modes. The upper panel is for the case of 9b = 0° 
and the lower for 9b = 90°. Other parameters used are 
the same as those in Fig. 1. 
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Fig. 4. Surface emission flux spectra obtained with the full 
computation and with the diffusion approximation. Pa- 
rameters used are the same as those in Fig. 1. The Planck 
function of the same effective temperature is shown for 
comparison. 

4. Results and disscusion 

For a given effective temperature, the computed temper- 
ature profiles of the neutron star atmosphere with and 
without the diffusion approximation are plotted in Fig. 1. 
A zero-field temperature profile is also shown for compar- 
ison. Since opacities are generally reduced in the presence 
of magnetic fields, a magnetized atmosphere is more trans- 
parent down to a denser layer. The temperature profiles 
of magnetized atmospheres show a plateau near p ~ 0.1 
g/cm 3 . These plateaus are the manifestation of the differ- 
ent depths at which the O-mode and X-mode photons de- 
couple with matter. In general, opacities are larger for the 
O-mode photons and therefore they decouple with mat- 
ter in a region outer in the atmosphere than the X-mode 
photons. The temperature profile does not differ much 
for different magnetic field orientations if a full compu- 
tation without the diffusion approximation is performed. 
For those with the diffusion approximation, the deviation 
from the full computation is larger for the case with a 
smaller 0b- 

The beaming pattern, that is, the intensity in different 
directions of propagation, is complicated. It is a combina- 
tion of the effects of limb darkening and the anisotropy 
due to the presence of strong magnetic fields. It also de- 
pends on the photon frequency and polarization modes 
In Fig. 2 we show the beaming pattern of the two polar- 
ization modes for a 10 17 -Hz photon. Two cases, 9b = 0° 
and 90°, are shown. To show the frequency dependence 
of these patterns, the intensity spectra in certain propa- 
gation directions are shown in Fig. 3. In most cases, the 
X-mode intensity is generally stronger than the O-mode. 
It is also obviously seen in Fig. 2 that the diffusion ap- 
proximation ansatz, P v = -£-{u\, + -fxFi), is not a good 
description for radiations in a magnetized atmosphere. 
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In Fig. 4 computed flux spectra are shown. At a given 
effective temperature, due to the atmospheric effect, that 
is, higher energy photons have a smaller opacity, all spec- 
tra are harder than the Planck function of the same ef- 
fective temperature. In the presence of magnetic fields, 
the spectra get somewhat softer, but still harder than 
the Planck function. We particularly note that the dif- 
ference between the spectra computed with and without 
the diffusion approximation is noticeable and is distin- 
guishable with the current observation accuracy. Those 
with the diffusion approximation are generally softer than 
those without. The effective temperature will be over- 
estimated when fitting the observed spectrum with the 
model computed with the diffussion approximation. We 
therefore urge readers to be cautious when applying the 
diffusion approximation. 

The results discussed in this paper are those for a local 
patch at the neutron star surface. To compare with obser- 
vations, even for phase-resolved spectra, a global spectrum 
is needed, which is a combination of emission from differ- 
ent patches all over the neutron star surface. Gravitational 
light bending and redshift should be also considered. Since 
contributions come from different part of the stellar sur- 
face, different magnetic field orientations relative to the 
local surface normal are inevitable. Although such a com- 
putation has yet to be performed and knowledge of the 
temperature distribution over the whole stellar surface is 
required, our results (Fig. 2 and Fig. 3) already indicate 
that the global spectrum will be polarized to a significant 
degree in the extraordinary mode. More quantitative and 
detailed results will be obtained with a more comprehen- 
sive computation in the future. 
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